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ABSTRACT 


The problem of narrow band interference while transmitting broad-band signals 
like Direct Sequence Spread Spectrum is a common source of problems in Electronic 
Warfare. This can occur either due to intentional Jamming or due to unavoidable 
signal sources present in the vicinity of the receiver. Lack of improper information 
on these narrow band interferers makes it difficult to cancel them. 

In this thesis the above problem is addressed by using an adaptive notch filtering 
technique. Before adopting such a technique other methods like the Least Square 
Estimator and the Maximum Likelihood Estimator were explored. However the Kwan 
and Martin adaptive notch filter structure was found both relevant and suitable for 
the problem of interest. The Kwan and Martin method has the difficulty of increasing 
hardware complexity with number of notches. This makes it difficult to implement 
in real time. A new algorithm was developed for the purpose of implementing the 
structure in real-time. This new algorithm offers the same performance at reduced 
hardware complexity. This algorithm was simulated and the results were presented. A 
hardware feasibility is discussed by proposing a simple structure based upon existing 


commercial signal processing chips. 
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I. INTRODUCTION 


A. MOTIVATION 

The problem of suppressing narrow-band noise from a broad-band spectrum is 
of considerable importance ın the areas of Electronic Warfare and Anti Submarine 
Warfare scenarios. In this thesis a workable solution is proposed for the problem cre- 
ated by narrow-band interference. This solution will come in the form of an adaptive 
digital notch filter. However before the details of the solution to the problem are 
discussed, a brief description of the problem is given 


1. Electronic Warfare(EW) 


A typical EW scenario is given in Fig 1.1. In this scenario a transmitter 
transmits information over a long distance. As an Electronic Counter Measure(ECM) 
this information is coded and transmitted as a Direct Sequence Spread Spectrum 
(DSSS) signal. Inherently this signal is a broad-band sıgnal. However at the EW 
receiver there is often narrow-band interference from such things as push to talk 
systems(PTS) that swamp out the received DSSS signal. Sometimes for reasons of 
signal security PTS operating frequencies are varying with time. Under conditions 
such as these the EW receiver cannot function effectively. For proper functioning of 
the EW receiver we must enhance the received signal by selectively and adaptively 
suppressing these narrow-band signals. 

2. Assumptions 

This thesis addresses the problem arising at the tactical data communica- 

tion link. In this EW scenario as dipicted in Fig 1.1 we are attempting to perform 


signal analysis on a DSSS signal using a EW reciever. This EW receiver is not the des- 
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ignated receiver and hence does not have the spreading code available. For this reason 
the narrow-band interference severely limits the ability to determine the charaterstics 
of the recieved signal such as carrier frequency, chip-frequency etc. 

A typical DSSS system is the Tactical Information Distribution Systems 
(JTIDS)[Ref. 7]. Data on such a link is generally in the form of digital messages. 
A typical L-band (960-1215 Mhz) JTIDS allows transfer of digital data between any 
properly equipped users within line-of-sight. The typical RF band-width required is 
around 10Mhz. 

Assuming a Superheterodyne EW-receiver, the band-width of the Interme- 
diate Frequency (IF) need be only 10Mhz. This enables us to digitize the signal at 
the base-band with a (20-50)Mhz Analog to Digital(A/D) converter. It is assumed 
that this signal 1s digitized and converted into a floating point number. In this thesis 
it is assumed that the floating point arithmetic 1s of sufficient precision that rounding 


errors are not a significant problem. 


B. PROBLEM DEFINITION 


The signal described in section (1.A.3) can be thought of as a broad-band signal 
with additive narrow-band signals occurring at arbitrary frequencies and times. These 
narrow-band signals can be modeled as either pure sinusoidal signals, as outputs of a 
narrow-band filter excited by white noise, or as any band-limited narrow-band signal. 


We shall assume that the signal yx is received by the receiver and is defined as 


3 
yk — wk d 5 ZEF (1.1) 


s=1 
where z}, the interferer, is either a pure sinusoid or a narrow-band signal whose 
frequency (i.e center frequency) can vary slowly with time. The term wy is the 


broad-band DSSS signal of interest and +, is White Gaussian Noise( WGN) with zero 


mean and variance 0? N(0, 0?) The goal is to remove only the narrow-band interferers 


3 


and obtain the output (wk + yk) for further processing. A typical spectrum for yx is 
given in Fig 1.2. 

In this thesis, the design for an adaptive digital filter that will achieve the goal 
of removing narrow-band interferers from a broad-band signal of interest is developed. 
This thesis is organized into five chapters. chapter I is the introduction. The chapter 
II concentrates on building necessary background on such things as the Least Square 
Estimator and Maximum-Likelihood Estimator. Chapter III briefly explains existing 
methods of adaptive filtering with emphasis on adaptive notch filtering techniques 
. The new algorithm [Ref. 26] developed for the purpose of addressing the current 
problem of interest is also explained. Chapter IV is totally devoted to modelling 
various signals required for testing. Simulation results are discussed and presented. 
Chapter V provides a brief survey of existing hardware components and a look at 


possible hardware structures and their limitations. Chapter VI is conclusions. 





| | L 001 


s 


00% 


Ä USC 


00E 


OSE 
uni32ads Teubrts [eordA3 wv 


II. ESTIMATION TECHNIQUES 


A. DIGITAL FILTERS 
Digital Filters can be characterized by their impulse response, their trans- 
fer function, or by a difference equation [Ref. 13] . A typical Infinite-Impulse- 


Response(IIR) filter represented in difference equation form is given by 
Tk = ayzk-1i മാമാ asz,-3 + biu, + b2uk-1 (2:5 
while an Finite Impulse Response(FIR) filter is given by 
zk = bur + b2Uk-1 + b3Uk-3 (2:2) 


where uz is the input to the filter and x, is the output of the filter. z and u, are 
the sampled values of the continuous functions z(t) and u(t). By choosing T, as the 


sampling time we notionally write the discrete signal as 
Te SE (2.3) 


Let X(N) be the spectrum of the signal Ty, where Q is the angular frequency 
component present in the actual signal. For analysis purposes consider normalized 
frequency given as w = 27,. With this definition we do not have to use the actual fre- 
quencies but only the normalized frequencies. Assuming that there is no aliasing while 
sampling the maximum normalized frequency is 0.5 cycle/sample or x rad/sample. 

1. Some Hardware Schemes 

Fig 2.1 represents a conventional Direct Form II representation of a digital 
filter [Ref. 18]. However for high-speed throughput it is often necessary to have a 
pipelined architecture. FIR filters are highly amenable for pipelining thus giving a 


high throughput rate for computations. 
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2. Modelling Hardware Multiplier and Adder 

Multiplication of two variables say bi and u, in hardware can be modelled 

as 
bi*u, = b; *u,D, = bur, D) (2.4) 
where x is the hardware multiplier and * is the ideal multiplier and D, is the propa- 
gation delay associated with the multiplication. Similarly a hardware adder is desig- 
nated as + while the ideal adder is given as + and the corresponding delay associated 
with the adder is D2. Then an addition of two variables u, and ux-) Can be written 


as 
urtur-ı = (ur + uk-1i) Do (2.5) 
Each of the previous devices can be Incorporated into a pipeline digital filter structure 
by following ıt with a synchronizing register. The clock period must then be less than 
D +t, +t, where t, is the register setup time and t, the register propagation delay. 

3. Pipelining an FIR filter 
In this section we wish to realize 

Tk = b * u, + b; * uy 1 + b3 * Uk? (2.6) 
It is acceptable to have an overall delay D , However it is desirable to minimize N: 


£y — D^ (b, * uy -F bo * uy 1 4 bs * uy 2) (2.7) 


First associate one D with each multiplier and form real multipliers by the notation 


Dx = *. Then we get 
Zp = D^- {b D *u, + bsD *uk — 1-+ bsD * urz} (2.8) 


Now form the real multiplier 


^ 


$, = D^-! (buy 4 b; suk — 1 4. b3*uj 2] (2.9) 


Š 


Now we use one additional D to associate with the adding of b2*u,_1 and bs*uk-2. 
(Note: It is important to choose the most delayed values first in order to minimize 
N) 

ž, = DM"? (Dbjžu, + D(božuko į + bažux-2)) (2.10) 


Now form the real adder 
£,— D^ [Db uy * (bo*us ai 58, 2)] (COSTI 


Now we use an additional delay D to associate with the final adder to get a real 


adder: 


fx = DN f Dbytur+(b2tun-ı+bstur-2) } (2.12) 
Finally we choose the system delay such that D = 27? and replace u4-2 = 
z-^u, and ui 1 — z^ lu,, thus we get 
jQ,2z 09 (2 biturt (27 btur tz batur) } (2315) 
Now factor z^! out of the expression to obtain 
f = z7 N+? {biturt (batutz  b3*u,)) (21) 


Note that N = 2 is the minimum possible value of N for a causal filter. Choose 


52 then 
m = (biku,+(b;*u,+z lbstu,)) (lo 
Equation (2.15) is very convenient for pipelining and is given one to one in Fig 2.2. 
The above procedure 1s general for any FIR filter. Pipelining an FIR filter 


results in maximizing the sampling or throughput frequency at the expense of a delay 


or latency in the availability of the output. 
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4. Pipelining IIR filters 


Pipelining an IIR filter is considerably more difficult than an FIR filter. 


This is mainly because the delay in the availability of the output makes it impossible 


to feed back output values witlı short delay as is required for the straightforward IIR 


implementation. As an example consider a simple 1st order IIR filter 


Th = 0] * Tk-1 T Uk 
As in the FIR case we introduce a delay D" as follows. 


^ 


£y, — D'(a,* zy 1 4 ux) 
Now we use one delay to form the real multiplier: 


tk = D~ (aD * zk- + Dux) 


DNajéxy-1 + Dux) 
We use one more D to form the real adder 
Tk = Ba * Tk-ı tur) 


Finally assume D= z’! and x,-ı — z^!z, then 


ടി) 


m oz taz, tu) 


For tą = zk, N must be zero. However N > 1 is required for a causal filter. 


(2.16) 


(2.17) 


Solutions to this problem exist. They use a higher order difference equation 


representation of the filter with equivalent characteristics, so that feedback loop delay 


greater than 1 can be tolerated. More details can be found in the literature [Ref. 14) 


regarding pipelining of IIR digital filters. 
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B. SIGNAL CHARACTERIZATION 
The signal y, of equation (1.1) can be characterized either in the time domain 
or frequency domain and the signal can be easily transformed from one domain to the 
other via a linear transformation. 
1. Time Domain Representation 
In the time domain representation, the signal yx is modelled as the output 
of a linear system which may be either an all pole or a pole-zero system [Ref. 13]. 
The input is assumed to be white noise but this input to the system is not accessible 
and is only conceptual in nature. Let the system under consideration be 
p q 
Yk = > 2 SE 22 bite; (2.23) 
where x is a white noise and yx is the output of the system. Equation (2.23) also 


can be represented as a Transfer Function(TF) (Ref. 17] 





E) Aa (2.24) 
Tk = | Aa | Vk (2.25) 
where 
A(z) = 1 — Dar (2.26) 
and 


q . 
Pi > bau (2.279) 
= 
we define the parameter vector as 
p = [21,42,---,a,5,b;, 62, , b] (2.28) 


The parameter p completely characterizes the signal 1; The above system defined 


by the equation (2.23) is also known as an Auto-Regressive-Moving-A verage( ARMA ) 
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model [Ref. 12],[Ref. 15]. It is conventional to denote a p poles and q zeros system 


as ARMA(p.q). 
2. Frequency Domain representation 


The signal yx of equation (2.23) can also be represented as a vector 


Y y = ly1 +++ yn] (2.29) 


which can be transformed into the frequency domain by the DFT relation 


) (N 
= | = (wya) } (2.30) 


n=0 
where W — e- and J = /—1. In general gj is a complex quantity. However we 


restrict our interest to the power spectrum a real quantity given as 


Sk = gk9k 231) 


The series SN 

SN = [sie sn] 232) 
is another form of representing the signal and ıs designated as the power spectrum of 
the signal z,. Even though the signal in this domain is very convenient to handle, 


computational complexity and other considerations limit its use for online application. 


Also in this representation phase information of the signal 1s lost. 


C. SPECTRUM ESTIMATION 
1. Least Square Estimator 


Consider the signal yx the output of a linear system excited by a white 
noise as given by the equation (2.23). Now our problem is to estimate the parameter 


vector p by obtaining successive measurements of yx. We shall define the vector 


Xi, = [Yk-1, tts Yk-ps Yks’? al (2.33) 
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and zx = yx so that the difference equation (2.23) can be rewritten as 
Zk = PLEK (2.34) 
an elegant solution [Ref. 2] for the above problem can be written as 
Pr = Pe-1 + Hexxes (2.35) 
where 
: 21 
H, = p xa) (2.36) 
and 
Ck = zk Dx (2.37) 


The block solution for the eguation (2.34) can be also written in the form 


Px = Hy p <=] (2.38) 


ı=1 
Equation (2.37) can also viewed as Ar, — B+, and this is represented in Fig 2.3. 
Computing the matrix H; via equation (2.36) is computationally inefficient 
and a recursive solution via the matrix inversion lema [Ref. 3] is preferred and is given 


by 


v 


(2.39) 


H- SH 
ത്‌ -| k-1Xk X, I1, | 


l+x7Hı-ı% 
Staticians refer to the H, matrix as the covariance matrix while optimization people 
refer to H, as the Hessian of the objective function, where the objective function is 
given by 

J, = > e? (2.40) 

i=] 

Appendix A gives a computer program based on this approach that was used for this 
thesis work to identify system parameters via Eguations (2.35), (2.36), (2.37), (2.39). 
The above method is known to work well as long as the measured value of the state 
of the filter is free from noise. In the event z, is not noise free, estimates are known 


to be biased. 
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2. Maximum Likelihood(ML) Estimator 
The ML estimator for the above problem is conventionally obtained by 
considering a log-likelihood function [Ref. 16] and minimizing it. However for getting 


a better physical understanding the approach given by Young[Ref. 2] will be followed. 


EE). Z 


Consider a system as 





and an auxiliary system referred to as the model 1s given as 


I B(z) 

Tk = | 777 2.42 

j el Ë "7 
and we define the function to be minimized as L — ez(p) where e, — y, — £,. 


To minimize this function the derivatives and the Hessian of the function L where 
L = (yx — tx)” need to be obtained. First take partial derivative of L with respect to 


each parameter. For example 


ðL . O OT; 
Ob, 2 2( Yk => TN Ob, ) a OC Ob, ) (2.43) 


The partial = can be obtained by differentiating equation(2.42) as 


Q B(z) ച O bi 02 : 
a”) Ib: (| 12) (2.44) 
1 
A | 2.45 
[x] ^ (25 


Similarly the other partial rn can be obtained as 





0 (BE) | (ASE) - | 
Ba, Lan ( [A(z)] 5 2 
= — Ad, Yk (2.47) 





since 2 (B(2)) =0 as B(z) has no terms containing a4. Using equations (2.45) and 


al 


(2.47) the partials are: 





el 
0 ET A(z) 








Oi, š B(z) 
da 7*7 E. 
Our k 1 

da, KS ir dh 


by defining 
t 
Pi = [a,, az, az, by, be] 
t * * * * * 
x — [pm Tk_2, Tk-3, Up» u; 1] 
t = = = = se ] 
Zk = |Yk-1>Yk-2>Yk-3>Uk>Uk-1 


the Hessian may be approximated by: 
L -1 
H, = p xat) 
1-1 
Using the above equations the parameter can be estimated via 


Pk = Pı-ı + uU Hixiey 


(2.48) 


(2.49) 


(2.50) 


(2.51) 


(2.52) 
(2.53) 


(2.54) 


(2.55) 


(2.56) 


where u is the step size in incrementing the parameter vector. A block schematic for 


the algorithm was given in Fig 2.4. A serious drawback of this method is the stability 


while incrementing the parameters. 


3. Other Methods 


There are many other methods such as the Instrumental Variable (IV) 


approach, the Approximate-Maximum-Likelihood (AML) Method and a combination 


method IV-AML available in the literature [Ref. 2]. 
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D. SELECTION OF ESTIMATION TECHNIQUE 

Before we select the appropriate method for the problem of ınterest, model-order 
problems must be considered. 

1. Model Order Problems 

In the Least-Square Estimator for the parameter vector p of the signal 

zk we need to assume the dimension of p or the order of the system. Since the 
system order is often not known a higher-order model estimate is assumed. Now we 
investigate higher-order model fits for a lower-order data. 


Consider the system equation (2.23) with parameter given as 
p! = [1.7, —1.53, 0.648, 1, 0.6] (2.57) 


For the above system a Random Binary Noise(RBN) was given as the input uz and 
the corresponding output z; is obtained. For the time series uy and z; an ARMA(6,1) 
model was fit using the program given in the appendix and the parameter was esti- 


mated as 
p' = [2.26533, —2.79577, 2.17249, —1.09877, 0.446598, —0.110151,0.998175] (2.58) 


Fig 2.5 gives the parametric spectrum for the above. 
For the same data an ARMA(4,1) model was used and generating a pa- 


rameter vector 


p' = [2.16988, —2.41286, 1.47436, —0.365412, 0.999317] (2.59) 


Fig 2.6 gives the parametric spectrum for the above values. 
Lastly for the same data an ARMA(8,1) model was used and generating a 
parameter vector 


p' — [2.276, —2.835,2.267, —1.282,0.703, —0.353,0.143, —0.0035,0.9999] (2.60) 
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Fig 2.7 gives the corresponding spectrum. 

It is interesting to compare these with the actual parametric spectrum given 
in Fig 2.8 For the sake of completion and comparison, the spectrum given in Fig 2.9 
of the output time series rk computed via a DFT program is also included. This 
simulation demonstrates that the existence of multiple solutions, hence an important 
step in the estimation procedure is to choose an appropriate model order. 


2. Choice of the Method 
Since the filter hardware must be implemented in real-time at the frequen- 
cies of 10 to 20 Mhz computational complexities must be kept to a minimum. These 
methods although providing good performance, are not well suited to hardware im- 


plementations. This motivates looking into new filtering methods. 
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Ml. METHODS for FILTERING 


A. FILTERING TECHNIQUES 

Filtering in a broad sense is selectively suppressing a portion of the spectrum 
of the given signal. In this chapter several filtering techniques are explored in the 
light of the narrow-band interference problem in order to identify applicable filtering 
techniques. 

1. Filtering via FFT 

Any given signal can be conveniently transformed into its power spectrum 

via equation 2.31. Assuming the spectrum of the desired filter to be defined by wi. 


The weighted output 1s given by: 
UK = SWK (3.1) 


and the corresponding filtered output y; is obtained by taking the inverse DFT of the 
signal vy. A simple block schematic representing this idea 1s given in Fig 3.1. Details 
of this approach are available in various references [Ref. 22]. 

Fast algorithms such as the FFT for computing the DFT make it possible 
to do the above process in real time by using dedicated hardware. Honevwell makes 
the HDSP66110 and HDSP66210 Digital Signal Processing chip pair which are ideally 
suited for these applications. This chip pair can perform a single complex multiply in 
40ns [Ref. 28]. However the problem of filtering adaptively [Ref. 10] (t.e varying wx 
according to a criterion) demands more computing power which can be obtained by 
augmenting the DSP chip pair with a processor like the R3000 which has a Reduced 
Instruction Set Computer( RISC) architecture with a high instruction execution rate 


of approximately 25 Million Instructions Per Second( MIPS) [Ref. 27]. 


23 





| 
| | 
| 
| n e. å Io — (—- Cx A | 
| X, | Uuonoung As | x x 
| " | | L |. ^OpUIM D TI 
"q + er Buiyybiam | war ` | 
| | m el |. A | 


L44 ela BuliƏ1l!3 








2. Recursive DFT 
Implementation of a higher order FIR filter using a recursive DFF ef. 
23],[Ref. 24] is also very convenient for eliminating narrowband interference. The 
basic concept is that the FIR filter is expressed as a product of two filter sections. 
One section is a filter with its zeros being equally spaced on the unit circle. This 
is achieved by a delay. The second section is a pole-producing section. Pole-zero 
cancellation results in the desired FIR filter. More dctails can be seen in the references 
[IRef. 23],[Ref. 24]. 
3. Adaptive Filtering 
Adaptive filters can be placed into four classes based upon the choice of 
the training sequence and the reference model for adaptation. Simon |Ref. 21|has 


classified the Adaptive Filters into the following four classes 


e Identification - Class I 

e Inverse Modelling - Class I] 

e Prediction - Class III 

e Interference Canceling - Class IV 


Fig [3.2 - 3.5] give the block schernatics of the various classes of adaptive 
filters. In our problem, the reference model is naturally a narrow-band bandpass filter 
since the interference signal is a narrow-band signal. The class of filtering that best 
suits our problem is a combination of Class III and Class IV type of filter [Fig 3.4 
and 3.5]. This logically points to an adaptive notch filter. Duc to apriori knowledge 
of the interference signal, only the parametric approach was. However the DFT- 
based techniques may also offer a good solution and should be the subject of future 


Investigations. 





p 
C — 
1indino uieisAiS 


1ndino wejsÁs 





ee Did 


Áejeq 


191114 
n 


e^] depy 


( 1| $8912 ) 6urjopow ƏsuəAu| 


ce 513 


< 
p 
301113 


oAlide py 





> 


(| SSe/9 ) uoneOSljnuəpi 


| 
| 





Indu} 
wejsÁs 


ınduj 
wejsÁs 


26 





5 101113 


e^ de 
indino wejsÁs Á CEER " Igu6IS @oSuele)əH 


Igu6IS A4BUuldd 


( Al S8810 ) B5uilləƏSueO ƏouəƏlJəjiəƏlui 


ve ലച 





Įgu6įs 
WUOPUBH 





103113 
e^ depy 





| indino uejsAÁg 9 


Z yndjno weysks 


( M4 $8919 ) uonotpeJd 


D 


B. ADAPTIVE NOTCH FILTERS 

Notch filters for removing multiple narrow-band interference can be categorized 
into four broad categories illustrated in Fig 4.6 through 3.9. The first two categories, 
Figures 3.6 and 3.7, are cascaded second order notches with each second-order section 
removing one freguency. The next two categories, 

Figures 3.8 and 3.9, are higher-order notches that eliminate multiple freguen- 
cies. In all of the categories, it is possible to use FIR filters (t.e all zero filters) which 
are easily pipelined and can be made truly linear phase. However, IIR notch filters 
require substantially fewer multipliers and adders than FIR notch filters. Thus IIR 
pipelining may become an important issue. In this thesis we limit our discussion only 
to the first two categories illustrated in Figures 3.6 and 3.7. 

1. Second-Order Cascaded Notch Filters 

The second-order notch filter is used in cascade and in-line with the signal 


as shown in Figure 3.6. The transfer function for such a notch filter is given by Kwan 


and Martin [Ref. 8] as: 


Hx(z) = 1- TA (3.2) 
ką k; (1 + z-1)(1 — z!) 
r dm | 2 1-— (2 = k; E Ka + (1 = ka. a 
പ 1- DR pa (3.4) 


uer eeu DU aU 

For arbitrary values of k, and kz , this is a symmetric notch filter with unity 
gain at DC and the Nyquist frequency. If kz is kept constant, then the 3db notch 
width is also kept constant. Thus k, may be adapted to remove one narrow-band 
signal. A cascade of such filters can be used to remove multiple narrow- band signals 


[Ref. 8]. Constants kı and kz are related to the pole radius r and the normalized pole 
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frequency 0, as 


ko = (1 = r°) (3.5) 
ky = 15-ന്‌? - 2rcos0, (3.6) 


It is important to bear in mind that the above transfer function has unity gain at 
i 
(peak = Zarcsin y — === (337) 
2/1 - 2 
2 
which is different from 0, [Ref. 8] 
2. Second-Order Cascaded Signal Canceler 
The cascaded second-order signal canceler approach shown in Figure 3.7 
has the advantage that the desired signal does not pass through the adaptive filter. 
Instead, the band-pass filter is used to detect the narrow-band signal which is then 
subtracted from the desired signal. A constant 3db bandwidth notch can be achieved 


by selecting a band-pass filter with the transfer function: 





Hel) = Rr (3.8) 
NG) 
= (3.9) 


where D(z) is the same denominator as Hy(2) in equation (3.4) The signal-canceler 
structure is also nice for adaptation because it is relatively easy to generate sensitivity 
functions which are related to the gradient of Hi,(z) with respect to the frequency 
parameter ki. Figure 3.10 shows the block diagram of an adaptive version of this 
filter. The sensitivity function s(n) is obtained by differentiating the error signal 
e(n) with respect to the parameter kı. This can be easily derived by noting that 
e(n) = x(n) — y(n) and 
de(n) Oy(n) 


ak, = 0 (3.10) 
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to get the above partial derivative we consider 


y(n) = Hu(z)e(n) (3.11) 


the above eguation (3.12) can be easily differentiated with respect to k; by using the 


eguation (2.47) to obtain: 








y(n) _ _ | N(z) 
Ok, DE 


by recognizing that Hy,(z) = D from equation (3.12) the above equation can be 


7 r (n) (3.13) 














written as 
Oy(n) _ DS 4 
and by defining 
Die 

HE DG) (3.15) 

we get the sensitivity function as 
Simi Le = H,(z) = Hop(z)Hss(z) (3.16) 
1 


The equation (3.16) is given as a block schematic in Fig 3.10. The parameter k, may 


then be adapted by the formula: 


e sia m (3.11) 





a. Computing ||s(n)||“ 
Ideally ||s(n)||^ can be calculated by the relation 


n 


ls) Ts (3.18) 


t=n— N 
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The above scheme computes the average of the past N samples. However a weighted 
average of the past samples with highest weight for the recent sample can be done 


quite easily by a first-order low pass filtering [s(n)]^ as follows: 
va Vd HO — A)s?(n) (3.19) 


However the above lowpass filtering can also be carried out by a second order filter 


such as: 


bil = 22°" +27?) 
pia 2rıkz! + 772 


| Is(n)]? (3.20) 


s E — 2cos(50)27* + 27?) 


A | Įs(n)] (3.21) 


Yet another way to estimate ||s(n)||“ i simply: 


ese 
Uc 2 [s' (n)] (3.22) 


the above form is defined as zero order forgetting. It should be noted that equation 


(3.21) was used by Kwan and Martin [Ref. 8]. 


C. MULTIPLE NARROWBAND SIGNAL SUPPRESSION 

The second-order implementations of section (1.) and (2.) offer considerable 
advantage both in hardware complexity and in adaptive performance for multiple 
narrow band signal suppression [Ref. 8],[Ref. 25),[Ref. 26]. 

1. The Kwan and Martin Filter 

In a recent paper by Kwan and Martin [Ref. 8], the problem of detecting 

and enhancing sinusoidal signals in the presence of noise is addressed with a cascade 
of IIR adaptive notch filters which are used to dune the sinusoids. Each of the 
sinusoids is eliminated by a bandpass filter whose output is an enhanced version of 


one of the sinusoids. Hence this remarkable structure can perform both tasks with a 
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single adaptive filter configuration which is shown to be highly robust and performs 
extremely well. 

The major disadvantage of the Kwan and Martin structure is that the 
number of biquad sections needed in the adaptive filter configuration ıs given by 
N(N+3)/2, where N ıs the number of sinusoids to be detected and removed. This 
becomes impractical in real-time situations with more than 4 sinusoids due to the 
geometric increase in the required hardware. 

a. Kwan and Martin Structure 

The Kwan and Martin structure consists of a cascade of IIR notch 
filters one stage of which ıs shown ın Figure 3.11. Each stage consists of a bandpass 
filter with zeros at DC and the Nyquist frequency and unity gain at its peak frequency 


w;. Such a filter would have the following z-domain transfer function 


Per: Ta 
H = — = —— 323 
di 2 1 —2r;cos@;z71 + poo ( ) 
where 
r; — pole radius of the 1-th section 
O, = 272 |W, 
w; = peak frequency of the i-th section 


ws= sampling frequency 


Kwan and Martin identify two different methods for adapting the 
filter. Most of their derivation is based on what they call the constant bandwidth filter 
in which the pole radius r; is a constant and only the frequency w; is adapted. An 
alternative approach which keeps a constant Q is also discussed in Kwan and Martin 
paper. In addition, Kwan and Martin select the adaptive quantity in such a Way that 


it 1s fairly easy to determine the notch frequency from the adaptive parameter. From 


35 


Figure 3.10, we see that the notch filter for each section is the difference between 1 


and the bandpass filter, hence: 


HUE Hj (z) (3.24) 


b. Calculation of the gradient 
The basic structure of the Kwan and Martin adaptive filter shown in 
Figure 3.12 is a cascade of N sections of the form of Figure 3.7. The overall transfer 


function 1s given by 


N 
T(z) = IIHN() (3.25) 


+=! 
N 


= II(!- H;) (3.26) 


Kwan and Martin choose as their objective function J(z) the square 


of the output of the final stage of the cascade: 


a 
Ts 
N 
— 

| 


(E, (z))° (3.27) 
[T(z)]°(A(2))” (3.28) 


Hence, the gradient of the objective function J(z) is given by 


8J(z) 
Ok; 


OT (z) 
Ok; 





= 2E\(z)X(z) (3.29) 


Thus in order to find the gradient with respect to the adaptive pa- 


rameters k; , we must take the partial derivative of T(z) with respect to each k; 


3 Hi (2) 
= TH (2) - an (3.30) 


ı=] 


122 
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From equation (3.24) we have 


OH}, (z) Ə[1 - H;,(z)] 





= — | T 33l 
Ok, Ok, (3.31) 
OH} (z) 
= ee 332 
ðk; (3-32) 


Using the rule for differentiation as given by the equations (2.45) and (2.47) and from 


equation (3.23) we have 





OH. (2) = j —2k,27? 
= Hip(2)H3,(2) (3.34) 


where D(z) is the denominator of Ee HÌ, is given by the equation (3.15). 
Substituting equation (3.34) into equation (3.30) we obtain the over 


all sensitivity for the jth kı parameter as 


OT (z) 


നിം TES (=) Hj, (z) H3,(z) | (3.35) 


=] 


12 
j-2 N 
= Ha n Hy ( e) Ai (2) H3 (2) (3.36) 


By recognizing that 


N x 
= | II i (989) 


= 


and by substituting equation (3.37) in the equation (3.36) we get 


sig 





- [Hue | (z) (3.38) 


Figure 3.12 shows the Kwan and Martin realization of the complete 
adaptive system. The difficulty in the Kwan and Martin approach is in the generating 
of the product of notch filters without the notch filter 7, as required ın equatıon (3.36). 


To generate thıs product for each section, would require N — 1 biquads per section 
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resulting in a total of N(N — 1) biquads just to generate the product. Kwan and 
Martin are able to reduce this by using the output of the bandpass filter as the input 
to their cascade via equation (3.37). Since this output already has (j — 1) of the 
required H}.(z) factors in it, only (N — 7) additional biquads are needed for a total of 
0.5(N? — N). Adding this to the N biquads required to realize the cascade of notch 
filters and the N biquads required to realize the H?(2) factors, yields a total number 
of biguads given as 0.5(N + 3)N. 
2. The New Structure 
Figure 3.13 shows the improved adaptive notch filter structure [Ref. 25],[Ref. 
26]. The key to the improvement is the recognition that the output Eı(z) = T(z)X(z) 
for the cascade of the notch filters can be written both as a product of the individual 
notch filter section transfer functions Hi (2) times the input and in terms of the input 


X(z) minus the outputs Y'(z) of the bandpass filters 


E) (z) TA) (3.39) 


N 
Take] X(z) (3.40) 


N 
X(z)- p v) (3.41) 


To get the product of Hi,(z) without the term ? — j as required in the 


equation (3.36), we may use equation (3.41) to simply add back the term Y?(z) 


N N 
IE) X (z) X(z) — 2,Y'(z) (3.42) 


EF EN) 
Ej (2) 4- Y?(z) (3.43) 


Figure 3.13 makes use of this fact to generate the gradient needed for the 


39 


പല. Did 





ss dq 
(Z) ya (Z) „N 

SS dq 
te 
ež 08 


mi p" 
SS Sr dq 
Is 3 6 


SINIONAS MON 


wa 





3 


(2) A 
N 


40 


adaptive process. From the Figure 3.13, we can see that the total number of biquads 
required is N for the cascade of notch filters plus 2N for the Hj,(z) H(z) required for 
adaptation, minus 1 at the last stage, since the last stage does not need the extra 


Hy {z)). Thus we have the number of biquads required in the new structure is (3N-1). 
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IV. MODELLING and SIMULATION 


A. MODELLING OF SIGNALS 
In order to test various algorithms and to evaluate their probable performance 
in the real environment, it was necessary to develop a meaningful simulation of the 


real situation. To achieve this, the following four testing categories were developed: 


(1) Sinusoidal signals with white gaussian noise 
(ii) Narrow Band Noise with white gaussian noise 
(iii) Bi-Phase Shift Keying (BPSK) sequence 
(iv) Frequency Shift Keying (FSK) sequence 
1. Sinusoidal Signals 
In order to generate sinusoidal signals with minimum computational burden 


placed at different frequencies 5 fs, a second order AR process 
zi = 2cos(0;)zi , — zi, (4.1) 


was used. Initial conditions are very important and they are chosen such that z_, = 0 
and r_2 = —sin(0,) giving a unit amplitude sinusoidal signal. The 0; value is between 
0 to 180 and n is the number of frequencies desired. The required signal y; needed 


to input into the adaptive algorithm is given as 


Vk = 2 Zk + % (4.2) 
t=1 
where yx is a white gaussian noise N(0,0?). 


2. Narrow Band Noise 


The narrow band noise signal is generated using the difference equation 
zi = 2rcos(;)zi_, riel + ul (4.3) 
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where 6, decides the placement of the noise in the spectrum and r controls the band- 
width of the noise. The u! is simply a uniformly distributed noise taken at different 
instants. The desired signal yx is obtained via 
izn 
Vk — 9 Zk 0 (4.4) 
i=) 
3. Bi-Phase Shift Keying Sequence 

The generation of BPSK signal has three distinct three parts: 

(1) Generation of Random Binary Sequence(RBS) is achieved by passing 
uniformly distributed noise through a hardlimiter (An important note is that the 
interval between the two consecutive bits of RBS is m 

(11) Another sequence of binary numbers is a spreading code or sequence. 
The specific sequence used in a given communications system 1s normally not available 
to anyone but to the designated receiver. (In this particular simulation we have 
generated the spreading sequence by passing uniformly distributed noise through a 
hardlimiter. The bit interval is m) 

(111) Phase encoding (t.e.) mapping the given binary signal which is the 
exclusive or of (1) and (11) as 0 or m at appropriate sample time. 

The output of the first hardlimiter is stored in an array z. Output of the 
second hardlimiter is stored in the array y This information is retrieved by a subtle 
use of the array index given as ? — kf, where ? is the index, k is the discrete sample 
number and f, is the bit rate of the intelligence. Similarly another index j is generated 


using J = kf. where f. is the chip frequency. Generation of the indices is the key 


thing in this simulation. The desired signal now is given by the equations 
wi = Acos(27 fok + à) | (4.5) 


$k = [z (i) @ y(g)]z (4.6) 
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p=n 


a, = > wi, (4.7) 
p=1 


where 2 and j are indices of the arrays as defined earlier. lt is an important point 
to note that the characteristic of the signal w4 are defined by the parameter p = 
(fo Je, Jo): | 

This signal is not really a simple BPSK signal but it has an additional fea- 
ture of spreading the spectrum by controlling the chip-frequency and carrier frequency 
and information rate. À block diagram of the scheme is given in Fig 4.1. 


The desired signal y; 1s given by: 
Yk = Qk + Pk (4.8) 


where o, is the set of narrowband BPSK signals placed at different places in the 
frequency spectrum and ( is the broad band BPSK signal generated for a specific p 
value. 
4. Frequency Shift Keying Sequence 

Generating this sequence needs a random binary intelligence signal. This 
was once again is achieved by passing a uniformly distributed noise through a hardlim- 
iter. The output of this hard limiter stored in an array z. An index z is chosen such 
that 7 = kf, where f, is the baud-rate of the information and k is discrete sample 


number. Now the desired signal is generated via 


s, — 2cos(0,)sy 1 — Sk-2 (4.9) 
O, comit) (4.10) 
ട്ട (4.11) 


where @ is the carrier frequency and ¢ is the depth of the frequency modulation. 
Initial conditions are very important and they are chosen such that s_ı = 0 and 


s-2 = —sin(0) giving an unit amplitude sinusoidal signal. 
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B. SIMULATION 

The adaptive digital filter algorithm as described in the earlier chapters is sim- 
ulated and tested using synthetic data. Simulation was carried out on a VAX 11/785 
computer using Fortran 77. Listing of the program used is enclosed in the appendix. 


The adaptive filter parameters of interest are 


a) Sharpness of the notch filter defined by pole position (r? = 1 — kz) 


b) Step size in the parameter update procedure (p) 


(a) 
(b) 
(c) Time constant of the fading filter (À) ( refer to equation 3.19) 
(d) Model order (N= # of 2nd order filters) 

(e) 


e) Order of the incoming signal or number of interferers (m) 


1. Kwan & Martin Algorithm 
In simulating this algorithm, the most important thing is the implemen- 
tation of the notch filter. Fig 4.2 gives the structure of the algorithm. Let z(n) be 
the input to the Adaptive Filter and e(n) the desired output. This desired output is 


obtained by passing the input r(n) via a cascade of N notch filters as 
e(n) (൧)... ൧(൭)] ര) (6 


Output at the intermediate jth section of the band-pass filter is designated as y(n) 
as shown in Fig 4.2 Then the sensitivity of the x? parameter, is obtained by passing 


this y’(n) through another cascade of (7 — 1) notch filters given as 
sn) = (EN HN HEHE) vm) (4.13) 


The filter transfer functions H} (z) and H?, were defined in the earlier chapter by 


equations (2.21) and (2.15). 
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2. The New Algorithm 
The primary difference between the New Algorithm and the Kwan and 
Martin algorithm is in the method of obtaining the sensitivity of the jth coefficient. 
This difference can be seen in the Fig 4.3. The output y’(n) at the jth notch Fig 4.3. 
filter is added to the over-all output e(n) and this summed output is passed through 


a cascade of only two filters as shown in Fig 4.3 and can be given as 
s(n) = |H},H3,| {e(n) + y’(n)} (4.14) 


3. Forgetting Filters 


For parameter incrementation via equation (3.17) we need to obtain ||s?(n) ||. 


This could be achieved by 


a) Zero order forgetting using equation (3.22) 
b) lst order forgetting using equation (3.19) 
c) 2nd order forgetting using equation (3.21) 


4. Stability 


The parameter incrementation given by the equation (3.17) can be recast 


by using either of the forgetting filters given above as 
$ t nom 
Ky = Kj —pu[v] s'(n)e(n) (4.15) 


It 1s very important to note that the above equation (4.15) has a close resemblance 
with the ML Estimator. In the MLE case the stability while incrementing the param- 
eter is an important factor. A similar problem exists in the current incrementation 
procedure but the problem 1s solved by a suitable choice of parameters and modify- 
ing the incrementation procedure by adding an additional factor to equation (4.15) 
as follows: 


| A 


Umaz 


Ki = xj — ue(n)s'(n) | tna 
min 


(4.16) 
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This modification also protects against possible overflow or underflow while computing 
[v] ^. In addition simulations indicate that in the case of zero-order forgetting, this 
factor is required in order to obtain convergence. When Vmin ANd Umar are zero and 
infinity respectively, equation (4.16) reduces to equation (4.15). 

In this incrementation procedure two forms of the denominator polynomial 


of the H,,(z) were considered: 
D(z)=1-(2-k2-kiJz* + (1— k2)2* (4.17) 


and 


D’'(z)=1-2«Kirz! + r?z”? (4.18) 


where Ki = cos(0'). 

While using the incrementation procedure, under transient conditions poles 
of the D'(z) cross the unit circle causing instability problems. This condition was 
averted by checking the pole position after incrementation and if unstable then in- 
crementation is modified. Checking this condition for D'(z) given by equation (4.17) 
calls for solving a quadratic equation at every parameter increment and examining 
the pole position. However this check is very simple for D'(z) given in equation 
(4.18), since we need only to check x! by maintaining |x{| < 1. Furthermore the 
value of r must be positive and less than unity for stability. The choice of the r ıs 
very important for proper fast convergence. Fig 4.4. shows the frequency response of 
band-pass filter for different values of r The 3dB band-width of the band pass filter 


Is related to r and is given [Ref. 29] by 





2 
Band - width = cos"! t | (4.19) 


Note that under limiting conditions, the band-width ıs essentially zero. In thıs simu- 


latıon ıt was assumed a value of r as 0.95. 
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A large number of simulations [Ref. 29] were carried out using sinusoidal 
inputs with and without noise for different m and N values. In these simulations 
the denominator polynomial used was given by equation (4.18). The results were 
tabulated in Table 3.1. From Table 3.1 it is seen that the choice of the step size is 
an important factor and step-size must be optimized for a specific number of sections 
N. The values of Vmin and Umar did not pose any problem while incrementing using 
1st order or 2nd order forgetting filters. It seems that a 1st order forgetting filter 1s 
more effective than either zero-order or a 2nd order filter. 

After fixing the values of y and A the algorithm was tested using Narrow 
Band Noise signal yx for its performance. This signal was used only to tune the value 
of r (i.e. sharpness of the notch filter). The following simulation results are based on 
parameter adaptation for the D'(z) given by equation (4.17). 

5. Response for Sinusoidal signal 

A sample simulation output due to sinusoidal signal was shown in Fig [4.5 
- 4.7]. The input signal is composed of 3 sinusoids with normalized frequencies As 
A fs and 12 f. and WGN with o = 1. The spectrum of this signal is shown in Fig 
4.5. After passing this signal through the adaptive filter we could see that interferers 
were removed and only the noise was left behind. This is clearly shown in Fig 4.6. 
The adaptation process is shown in Fig 4.7. | 


6. Response for Narrow Band Noise 


In this testing category sinusoids were replaced by narrow-band signals. 
These signals were generated via equation (4.4). Superimposed on this signal, WGN 
with o = 1 was added. A typical spectrum of this signal is shown in Fig 4.8. This 
signal was passed through the adaptive filter and these narrow band interferers are 
notched out and only WGN is left out as shown in the Fig 4.9. Fig 4.10 shows the 


corresponding parameter convergence. 
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7. Tracking FSK Signal 
Transient behavior of the adaptive filter was studied by applying an FSK 
signal (generated via equation 4.11). In the case of an FSk signal, signal spectrum 
constantly changes. This property is useful for testing the tracking ability of the 
adaptive filter. In fact a WGN with o = 1 was also added to the signal. An adaptive 
filter then was used to demodulate the signal, by tracking the spectrum. This is 
shown in Fig 4.13. 
8. Suppressing BPSK Interference 
Having seen the performance of the adaptive filter under various conditions, 
it is only needed to test under an additional constraint of broad band noise. In this 
case, the broad-band signal was swamped by narrow-band interferers and WGN. A 
typical broad-band signal is generated by using a BPSK signal generator via equation 
(4.8). In this we have used carrier at 0.25 cycle/sample and the chip frequency at 0.125 
cycles/sample. The narrow-band interference signals are also generated by using the 
same BPSK signal generator, but with different parameters. Interferers are chosen 


such that the chip frequency is fixed at state cycles/sample, while carriers were chosen 


at 2 n fs, and ef, cycles/sample. To these ınterferers and broad-band signal, 
a WGN with e 2 1 was added. The composite signal spectrum is shown in Fig 4.14, 
indicating the three interferers, broad-band signal, and WGN. This signal was passed 
through the adaptive filter. Output of the filter is shown in Fig 4.15. In this figure 
it is clear that the interferers were removed by the filter and only the desired signal 
was left behind. 


9. Model order mismatch 


The model order of the algorithm was fixed at 3 while BPSK signals were 


generated with 2 interferers i.e. m = 2 and N = 3. Fig 3.16 shows the parameter con- 
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vergence under these conditions and Fig 3.19 shows for m — 1 and N = 3 conditions. 
We could notice that free extra parameter(s) wandering due to the mismatch. This 
has the effect of notching out the desired signal. This can be solved by deliberately 


introducing a known BPSK signal at fixed place. 
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Summary of Tests on the New Algorıthm 














TABLE 4.1 
Case | of | Sines | es | Umar | Noise 

T [TE [k [005 for [005 [m [0.05 — 
aa ജതി മങ്ക a oes He 
Hb |i |5 [005 20 [ne [605 [jump [0.05 [30 [005735 — 
[3b [5 5 [010 [10 [no 0.017 | 1410 | 0.063 | 150 [0.07 | 500 
(bs [5 |s 7040 [10 [yes 0.017 | 3500 | 0.063 [175 | 0.07 | 250 
He [5 — [r foto [10 [so [00171087 [0.063 [450 [0.07 | jump |. 
Den [5 [7] 0.10 [10 [ye [0.017 | 700 [0.063 | 167 | 0.07 | jump 
[3a [10 [7 [002150 [no [002 [166 [0.075 | 500 [0.07 | 1500 


In all cases with noise, parameters converge to exact value when number of 





notches is greater or equal to number of sinusoids. Iterations listed in table 


represent convergence to three decimal places. 


e In all cases with noise, parameter oscillate around the exact value. The number 
of iterations listed in the table is the number of iterations required to establish 


this oscillatory pattern. 


In some cases with more interfering sine waves than notches will jump at random 


between sine waves. This pattern is indicated in the table by the word jump. 


Values of Umin and Umar in the table refer only to zero-order forgetting which 
generally will not converge without limits on v. All other types of forgetting 


were run without limits on v 
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V. HARDWARE IMPLEMENTATION 


Although hardware design is beyond the scope of this thesis, in this chapter we 
shall briefly look at possible hardware configurations with a view toward feasibility 


of the new algorithm. 


A. HARDWARE FEASIBILITY 
In recent years many dedicated digital signal processing chips have become 
commercially available. Table 5.1 gives characteristics of some typical chips that are 


currently available. 


TABLE 5.1 









Feature Commercial Make 


and ത്‌്‌ 
[Oydefmuk 8 — [52 — s > j 
ease oe aH 


Data on MIPS R3010, Weitek 3364, and TI 8847 was obtained from the computer 








architecture book by Hennesey and Patetrson [Ref. 30] while the data for AT&T 
DSP32c is obtained from DSP32c data manual [Ref. 31]. Fig 5.1. gives a schematic 
of the 2nd order IIR filter hardware scheme. It has a coefficient memory which can 
be set by an external device. Close examination reveals that Multiplications A can 
be performed simultaneously while Afultiplication B and other additions are to be in 
sequence. Fig 5.2. is the basic building block for the adaptive filter and is identical for 
all the sections. This is offers an advantage in hardware complexity over the Kwan & 


Martin approach. Filters 1 and 2 have transfer functions of Hip(z) while filter 3 has a 


CK hotlok2 dela : 
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transfer function H,,(z). Throughput rate is primarily determined by the computing 
time of the 2nd order IIR filter. The overall block diagram of the adaptive filter for a 
3-Notch cancellation was given in Fig 5.3. This architecture remains same either for 
a Floating point or Fixed Point arithmetic. 
1. Time Budget 

In the architecture of Fig 5.1, 5.2, and 5.3 the most vital elements are the 
IIR filter, Controller and the Forgetting filter. The IIR fllter corresponds to H,,(:). 
The controller has to update the parameter via equation (3.17) which calls for 2 
multiplications, 1 addition, and 1 division. Similarly the Forgetting Filter has to 
implement equation (3.19) calling for 2 multiplications and 1 addition. Results are 


tabulated as given below 


TABLE 5.2 


Element + of Effective | 4 of Effective | # of 
u | kt | Discs 
RAR —— — [3 [m 
Controller? [1 E —— 


Fogetting 
Filter 


The maximum computing time is at the IIR filter. Using the AT&T DSP32c proces- 






sor, it takes 5x2 = 10 cylces (te 200ns) for implementation the desired IIR filter. This 
corresponds to a Throughput rate of 5M hz. This is the best that could be achieved 
with the existing commercial floating point processors. However processors like the 
HDSP66110 of Honey well make claim an ALU speed of 10ns with block floating point 
operations can conveniently implement at 10M hz throughput rate, without pipelin- 
ing the IIR filter. With pipelining the throughput rate would increase dramatically, 


but will never exceed the maximum throughput rate of 100M Az. 
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2. VLSI Approach 


Taking advantage of the fact that the entire hardware can be generated by 
a simple repetitive use of the building block, strong consideration should be given 
to designing a VLSI chip for Fig 5.2 rather than implementing the hardware via the 
existing commercial chips. The main rationale behind this idea is that reliability 


of the system is very important in EW equipment. Also VLSI has the potential of 


ultimate-low cost. 
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VI. CONCLUSIONS 


In this thesis the problem of suppressing narrow-band interference was ad- 
dressed. The problem specific to the Electronic warfare scenario was kept ın mind 
while solving the problem. The new algorithm [Ref. 26] was derived and simulated. 
This new algorithm was tested against various signals and signal conditions. The 
results are highly encouraging. 

Some aspects of pipelining were discussed. Reduced hardware complexity of the 
New Algorithm ıs an important advantage over the earlier algorithm by Kwan and 
Martin [Ref. 8]. A possible hardware scheme for implementing the new algorithm 
was discussed. It was observed that with the existing commercially available floating 
point processors a throughput rate of 5Mhz is achievable while using processors like 
HDSP66110 with an ALU speed of 10ns it is possible to achieve a throughput rate 
of 10Mhz. 


A. FUTURE WORK 


Future directions of work are 


1) VLSI design of the system with an architecture 
as indicated earlier or a similar architecture. 

11) Pipelining of 2nd Order IIR filter suited for 
this application 


111) Design using recursive DFT 


The above areas of work are the logical extensions of this work. However a radically 


different approach for the same problem using adaptive filtering in frequency domain 


[Ref. 10] should also be considered. 
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APPENDIX A 


dimension faray(5000) 

dimension uaray(5000),yaray(5000) 
dimension fhat(5000) ,ph(5000) 

dimension a(90) ,b(90) 

open (unit=9, file='spk.dat',status='new') 


This programme generates input output 
sequence by exciting a linear system 
defined by the numerator polynomial 

B(z) and denominator polynomial A(z) 
this data is stored in uarray and yarray 
Subsequently an ARMA (pole,zero) is used 
to fit this data. 


n=1024 
ix=1 
yk=rand(ix) 
ix=0 

kk=8 


generate sequence uk and yk 


ip=3 

iz=2 

a(1)=1.7 
a(2)=-1.53 
a(3)=0.648 
b(1)=1.0 
b(2)=0.6 

do 10 k=1,n 

call rbn(uk,ix,k) 
call system(a,b,uk,yk,ip,iz,k) 
uaray(k)=uk 
yaray(k)=yk 
continue 


save the sequence in the arrays uaray and yaray 
call spktrm(yaray,faray) 

spectrum of the output sequence yaray is computed 
using FFT programme and stored in faray for 


comparision 


ip=4 
iz=1 


Fit an ARMA(ip,iz) model for the given data 
this is an ordinary least squares algorithm 
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print *, "input-ip-12-->® 

read *,ip,iz 

do 20 k=1,n 

uk=uaray (k) 

yk=yaray (k) 

call arma(a,b,uk,yk,ip,iz,kk,k) 
print *, 'num->', (b(1),1=1,12) 
print *,'den->', (a(i),i=1,ip) 
continue 

print *,'num->', (b(i),i=1,12) 
print *,'den->',(a(i),i=1,ip) 


Sampling of the Z-transform. given 
coefficrents a .... and b ..... 
corresponding H(z) = B(z)/A(z) 

is evaluated on the unit circle 

for obtaing the parametric spectrum 


call zsmple(a,b,fhat,ph,ip,i2) 


do 30 i=1,n/2 

print *,i,'-->',fhat(i),faray(i) 
write (9,*) i,fhat(i),faray(1i) 
continue 

stop 

end 


subroutine zsmple(a,b,r,ph,ip,iz) 
dimension a(90),b(90) 
dimension r(2048),ph(2048) 
complex z,ui,fn,fd,omega,delw, spec 
pi=atan(1.0) 

pi=4.0*pi 

ui=(0.0,1.0) 
delw=(pi/1024.0) *ui*2.0 
omega=(0.0,0.0) 

do 100 k=1,1024 
z=exp(-omega) 

fd=(0.0,0.0) 

do 102121, 1P 
fd=fd+a(i)*(z**(-i)) 
fd=1.0-fd 

fn=(0.0,0.0) 

do 20 i=l,iz 
fn=fnt+b(1)*(z**(-1) ) 
spec=fn/fd 

r(k)=abs (spec) 

specr=real (spec) 
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speci=aimag(spec) 
ratio=speci/specr 
ph(k)=atand(ratio) 
omega=omega+delw 

print *,'--z-om-»',z,omega 
continue 

return 

end 


function atand(x) 
pisatan (1.0) 
pi=4.0*pi 
xrad=atan(x) 
atand=xrad*(180.0/pi) 
return 

end 


subroutine rbn(b,ix,k) 
z=rand(ix)-0.5 
if(z.ge.0.0) b=1.0 
if(z.1t.0.0) b=-1.0 
return 

end 


subroutine system(a,b,uk,yk,ip,iz,k) 


dimension 
& zkar(90) ,zkma (90), 
& a(90),b(90) 


if (k.ne.1) go to 250 
do 230 i-1,ip 
zkar(1)=0.0 

do 240 i=1,iz 
zkma(i)=0.0 
continue 

do 30 i=ip,2,-1 
zkar(i)-zkar(i-1) 
zkar(1)=yk 

do 31 1=iz,2,-1 
zkma(1)=zkma(i-1) 
zkma(1)=uk 


yk=0.0 

do 200 i=1l,ip 

yk=yk+a(1)*zkar(1) 

So 2191-1712 

yksyktzkma (i) *b (1) 
return 
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end 


subroutine delay(uk,yk,idelay,k) 
dimension d(500) 

if (k.ne.1) goto 10 
do 20 i=1,500 
a(1)=0.0 

continue 

do 30 i=500,2,-1 
d(1)-d(1-1) 

a(1)=uk 

yk=d (idelay+1) 
return 

end 


subroutine spktrm(taray,faray) 
complex u(5000) 

dimension taray(5000),faray(5000) 
m-10 

n=1024 

rn=n 

do 20 i=1,n 

u(i)=taray(i) 

continue 

call pfft(u,m,n) 

do 30 i=1,n/2 
faray(i)=real(u(i)) 

continue 

return 

end 


K k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k K 
* 

Subroutine for computing DFT of 

an array 'a' is complex and a pair 

of numbers are to be specified 


for each point * 

m is the 2 power index * 

say m=10 then number of points * 

are 1024 * 

after computation fft is kept * 

in the same complex array 'a' * 
* 


K k k k k kk k k k k k k k k k k k k k k k k k k k k k k k k k kk 


subroutine pfft(a,m,n) 
complex a(5000),u,w,t 
complex ui,ur 
pi=3.141592653589793 
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n-2**m 

nv2-n/2 

nml=n-1 

j=1 

do 10 i=1,nml 
if(i.ge.j) goto 20 
യി) 

a(j)=a(i) 

a(i)=t 

k=nv2 

if(X.ge.j) goto 10 

J=) K 

k=k/2 

goto 30 

j=j+k 

do 40 1=1,m 

le=2**1 

lel=le/2 

u=(1.0,0.0) 
w=cmplx(cos(pi/lel) ,-sin(pi/lel)) 
do 40 j=1,1e1 

do 50 i=j),n,le 

ip=i+lel 

t=a(ip) tu 

a(ip)=a(1)-t 
a(i)=a(1)+t 

u=u*w 

rn=n/2 

ui-(0.0,1.0) 
ur=(1.0,0.0) 

do 60 1=1,n/2 
a(i)=a(i)/rn 
a(i+512)=a(i+512)/rn 
areal=real(a(i)) 
aimgn=aimag(a(i)) 
amagn=abs (a(i)) 
aphase=atand (aimgn/areal) 
a(i)=ur*amagn+ui*aphase 
areal=real(a(i+512) ) 
aimgn-aimag(a(i-*512)) 
amagn=abs (a(i+512)) 
aphase=atand(aimgn/areal) 
a(i+512)=ur*amagntui*aphase 
continue 

return 

end 


$9$SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 


subroutine arma(a,b,uk,yk,ip,iz,kk,k) 
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dimension 
& zkar(90) ,zkma(90), 
& Zk(90),a(90),b(90),c(90) 
if (k.ne.1) go to 250 
n=ip+iz 
ipp=ip+1 
do 230 i=1,ip 
a(i)=0.1 
do 240 i=1,iz 
b(i)=0.1 
continue 
do 30 i-ip,2,-1 
zkar(i)=zkar(i-1) 
zkar(1)=ykold 
ykold-yk 
do 31 i=iz,2,-1 
zkma (i)-zkma (1-1) 
zkma (1) =uk 
do 32 i=1,ip 
c(i)=a(i) 
zk(i)=zkar(i) 
do 33 i=1,iz 
c(itip)=b(i) 
zk(itip)=zkma(i) 
call syseq(c,yk,zk,n,kk,k) 
do 34 i-1,ip 
a(i)=c(i) 
do 35 i=1,iz 
b(i)=c(i+ip) 

return 
end 


subroutine syseq(a,yk,zk,n,kk,k) 


EZEXZZZZXZkRERZEEEREREEREEEEREZEERZRERERREREEEE: 
solves system of equations with 

# of equations more than unknowns 
using linear reggression. 


equation 1S yk=zk * a 
where 'a! is n vector to be estimated. 


EXZEXZEXEZEREERSERERSEEREREEREREEREREERERZEREEREEE: 


dimension 
& del (90) ,phat (90,90) ,zk(90) ,y(90) ,a(90) 
if (k.ne.1) go to 250 
alpha=1.0 
continue 
call inv(phat,zk,alpha,n,kk,k) 
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ykhat=0.0 
do 40 i=1,n 
ykhat=ykhat+zk(i)*a(1) 
er=yk-ykhat 
do 60 i=1,n 
y (1) =er*zk(1) 
do 80 i=1,n 
del(i)=0.0 
do 80 j=1,n 
del(i)-phat(i,j)*y(j)*de1l(i) 
if (mod (k, kk) .ne.0) goto 70 
print *,'--->',k 
print 100, ((phat(i,j),i=1,n),j=1,n) 
print 130, (del(i),i=1,n) 
print 160,er 
print 170, (zk(j),j=1,n) 
continue 
do 50 i=1,n 
a(i)=a(i)+del(i) 
format (2x, 'phat',2x,/,4e16.6) 
format (2x,'2zk',2x,/,4e16.6) 
format (2x,'del',2x,/,4e1l6.6) 
format (2x, '!ek1'!,2x,el6.6,/) 
Format(2x,'zk',2x,/,4e16.6) 
return 
end 


subroutine inv(phat,zk,alpha,n,kk,k) 

dimension phat (90,90) ,zk(90),q(90),r(90) ,delp(90,90) 

if(k.ne.1) goto 5 

do 6 i=1,n 

do 6 j=1,n 

if(i.eq.j) phat(i,j)=1.0 

if(i.ne.j) phat(i,j)=0.0 
continue 
continue 

do 90 i=1,n 

do 90 j=1,n 

phat(i,j)=phat(i,j)/alpha 
continue 

do 10 i=1,n 

q(i)=0.0 

do 10 j=1,n 
alı)-phat(i,j)*zk(j)+q(i) 
sum=1.0 

do 20 i-1,n 
sum=sum+zk (1) *q(i) 

do 30 j=1,n 

r(j)=0.0 

do 30 i=1,n 
E(j)=phat(1,))*zk(1)+r(j) 
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do 40 j=1,n 
do 40 i=1,n 
40  delp(i,j3)-a(i)*r(j) 
do 50 j=1,n 
do 50 i=1,n 
50 phat(i,j)-phat(i,j)-(delp(i,j)/sum) 
return 
end 
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APPENDIX B 


simulation of michal paper IEEE 


dimension ak1(10) ,ak2(10),s(10) ,xout (10) ,theta (10) 
dimension aks(10) ,wp(10) ,ak3(10) 

dimension array (4,4000) 

open( unit=9,file='martin.dat',status='new') 
open( unit=8,file='mart.dat',status='new') 
print *,'-input-SNR-in-dB--->' 

read *,snr 

if (snr.ne.0) goto 200 

var=0.5 

goto 210 

var=1.0/(2.0*(10.0**(snr/10.0))) 

continue 

print *,'--variance-->', var 

n=1 

print *,'-input-nh--»' 

read *,n 

print *, '-input-#-of-waves-->' 

read *,nd 

print *,'-input-angles-->' 

read *, (theta(i),i=1,n) 

r=0.9 

pi-4.0*atan(1.0) 

rad-pi/180.0 

amu=0.0 

do 50 kl=1,50 

amu=0.015 

print *,'-input-step-size-- 0.015-->', amu 
read *,amu 


print *,'--angles-->',(theta(i),i=1,n) 

print *,'--input--1-for-step-change--else-0--»' 
read *,ic 

if (ic.eq.1) print *,'--input-step-in-degrees--»?' 


if (ic.eq.1) read *,step 

do 40 i-1,n 

ak2(i)-(1-r*r) 

ux-theta(i)*rad 

value=cos (ux) | 
akl(i)-di*sqrt(1.04r*r-2.0*r*value) 
akl(i)=sqrt(1.0+r*r) 
aks(1)=sqrt(1.0+r*r-2.0*r*value) 

continue 

amu=amu+0.0001 

do 10 k=1,6500 

call system(yk,wk,theta,step,var,ic,nd,k) 
call filters(ak1,ak2,aks,s,en,xout,yk,n,k) 
call increment (akl,ak2,aks,s,en,amu,avg,n,k) 
do 30 i=1,n 
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XX-2.0*sqrt(1-0.5*ak2(i)) 
wp(i)22.0*asin(ak1l(i)/xx) 

wp(i)-wp(i)/rad 

wp(i)-cos(wp(i)) 

continue 

do 60 i=1,n 

ak3(i)=abs(ak1(i)) 

continue 

kk=mod (k, 100) 

11=mod(k,10) 

if(kk.eq.0) print *,'--ak1-->',(ak3(j),j=1,n) 
if(ll.eq.0) write(9,*) k, (ak3(j),j=1,n) 
if(kk.eq.0) print *,'--aks-->', (aks(j) ,j=1,n) 
if(kk.eg.0) print *,k,'-avg->',avg 

print *,(wp(i),i-1,n) 

write(9,*) K,(Wwp(i),i=1,m) 
if(k.gt.1500.and.k.1t.2500) write(8,*) K,yk,en 
if(k.gt.5000.and.k.1t.6000) write(8,*) k,yk,en 
continue 

print *, kl,amu,avg 

write(9,*) kl,amu,avg 

continue 

close(9) 

stop 

end 


subroutine increment (akl,ak2,aks,s,en,amu,avg,n,k) 
dimension ak1(10) ,ak2 (10) ,aks(10),s(10) ,ss (10) 
dimension dec(10) ,flag(10) ,fading (10) 

if (k.ne.1) goto 30 

amu-0.001 
avg-0.0 
fading(1)=0. 
fading(2)=0. 
fading(3)=0. 
fading(4)=0. 
continue 

do 10 i=1,n 
ss(i)=fading(i)*ss(i)+s(i)*s(i)*(1.0-fading(i)) 
continue 

do 31 i=1,n 

print *,'-en-s-ss->',en,s(i),ss(i) 

sum=0.0 

do 20 i=1,n 

dec(i)=amu*en*s(i)/(ss(i)+0.0001) 
ak1(i)=ak1(i)-dec(i) 

sum=sum+ (abs (aks(i)-ak1(i)))/(aks(i)) 

call stability(ak1,ak2,flag,n) 
ak1(i)=ak1(i)+flag(i)*dec(i) 

continue 

realk=k 


0 00 0 
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aVg-avg*((realk-1.0)/realk)-*sum/realk 
return Ä 
end 


subroutine filters(ak1,ak2,aks,s,en,xout,xin,n,ki) 
dimension w(10,3),a(10,2),b(10,2),e(11),xbp(10) 
dimension qg(10,3),gbp(10),u(10) 

dimension q(10,3),s(10),xout(10),c(10) 

dimension aki(10),ak2(10),aks(10),as(10) 


1f (kl.ne.1) goto 200 

open( unit=9, file='martin.dat',status='new') 
continue 

i=n 

do 40 i=1,n 
a(i,l)=-(akl(i)*akl(i)+tak2(i)-2.0) 
as(i)=-(aks(i)*aks(i)+tak2(i)-2.0) 
a(i,2)=-(1.0-ak2(i)) 

b(i,1)=-0.5*ak2(i) 

b(i,2)=0.5*ak2(i) 

c(i)=-akl(i)*ak2(i) 

continue 

do 41 i=1,n 

print *,'-k1-k2-»',ak1i(i),ak2(i) 

pslnt *,'-a-b-c-5',aí(1,1),a(i,2),b(i,1),b(i,2),c(i) 


e(1)=xin 

i=n 

do 10 j=1,n 

w(i,3)=w(i,2) 

w(i,2)=w(i,1) 
w(i,1)=a(i,1)*w(i,2)+a(i,2) *w(i,3)+e(j) 
xbp(j)=b(i,1)*w(i,1)+b(i,2) *w(i,3) 
e(j+1)=xbp(j)+e(j) 

171-1 

continue 

൩-൦ (൩-1) 


5-൩ 

do 20 k=1,n 

u(k)=en-xbp(k) 
g(1,3)=g(i,2) 
g(1,2)=g(1,1) 
G (ny l)=a(1,1)#g(1,2)+*a(i,2)*g(i,3)+ü(k) 
gbp(k)=b(i,1)*g(i,1)+b(i,2)*g(i,3) 

q(k,3)=q(k,2) 

q(k,2)=q(k,1) 

q(k,1)=a(1,1)*q(k,2)+a(i,2) *q(k,3) +gbp(k) 
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Q0000 
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100 


10 


20 
10 


30 


s(i)=c(i)*q(K,2) 

i=i-1 
continue 
kk=mod (ki,10) 

if(kk.eq.0) write(9,*) ki,en,a(1,1),as(1),a(2,1),as(2) 
lf(kk.eq.0) print *, Kl,(a(1,1),1l=1,n) 

if(kk.eg.0) print *, ki, (as(i),i=1,n) 

if(ki.ge.1800) write(9,*) ki,e(1),en,xbp(1),a(1,1),as(1) 
if(ki.ge.1800) print *, ki,e(1),en,xbp(1),a(1,1),as(1) 
return 
end 


subroutine system(yk,wk,theta,step,var,ic,n,k) 
dimension theta(10),yout(10),yz(10) 
if (k.ne.1) goto 100 

continue 

call bpf(theta,yz,n,k) 

call waves(theta,yout,n,k) 

call file(yk,theta,step,ic,n,k) 
sum-0.0 

wk=0.0 

do 10 i=1,n 

sum=sum+yout (i) 

wk=wk+yz (i) 

continue 

call gnoise(g,var,k) 

yk=sum+g 

yk=wk+g 

yk=yk+g 

return 

end 


subroutine waves(theta,yout,n,k) 
dimension zkar(10,3),theta(10) ,yout(10) ,a1(10) 
if(k.ne.1) goto 10 
pi=(atan(1.0))*4.0 

rad=pi/180.0 

do 20 i=1,n 

zkar(i,1)=0.0 
zkar(i,2)=-sin(theta(i)*rad) 
al(i)=2.0*cos (theta (1) *rad) 
continue 

continue 

do 30 i=1,n 

zkar(i,3)=zkar(i,2) 
zkar(i,2)=zkar(i,1) 
zkar(i,1)=a1(i)*zkar(i,2)-zkar(i,3) 
yout(i)=zkar(i,1) 

continue 
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10 


20 


20 


10 


30 


return 
end 


subroutine gnoise(g,s,k) 
if(k.ne.1) goto 10 
ix-1 

yfl-rand(ix) 

ix=0 

continue 

sum=0.0 

do 20 i=1,12 
yfl=rand(ix) 
sum-sum-tyfl 
g-(sum-6.0)*s 
return 

end 


subroutine stability(akl,ak2,flag,n) 
dimension ak1(10),ak2(10),flag(10),a(10,2) 
do 30 is1,n 

flag(i)=0.0 
a(i,l)=-(akl(i)*akl(i)+tak2(i)-2.0) 
a(i,2)=-(1.0-ak2(i)) 
des=a(i,1)*a(i,1)+4.0*a(1,2) 

if (des.1t.0.0) goto 10 
des=0.5*sgrt (des) 
root1=0.5*a(i,1)+des 

rl=abs(root1) 

root2=0.5*a(i,1)-des 

r2=abs(root2) 

if (r1.ge.1.0) goto 20 

if (r2.ge.1.0) goto 20 

print *,'-real-roots--»',rootl,root2 
goto 30 

flag(i)=1.0 

goto 30 

des=-des 

y=0.5*sqrt (des) 

x=0.5*a(i,1) 

xx=abs (x) 

if (xx.le.1.0e-6) x=0.0000001 
angle-(atan(y/x))*57.3 
radius=sqrt(x*x+y*y) 


print *,'--complex--roots--»',radius,angle 
if(radius.ge.1.0) goto 20 

continue 

return 

end 
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subroutine bpf(theta,yout,n,k) 
dimension zkar(10,3),a(10,3),zkma(10,4),b(10,3) 
dimension theta(10),g(10),yin(10),yout(10),uk(10) 
if (k.ne.1) go to 10 
pi=4.0*atan(1.0) 
rad=pi/180.0 
do 15 i=1,n 
zkar(i,1)=0.0 
zkar(i,2)=0.0 
r=0.985 
print *,'--input--r-->',r 
read *,r 
do 25 i=1,n 
a(i,1)=1 
a(i,3)=r*r 
g(i)=0.5*(1-r*r) 
b(1,1)=9(1) 
b(1,3)=-9(1) 
a(i,2)=2*r*cos(theta(i)*rad) 
b(i,2)=0.0 
continue 
save=var 
do 35 i=1,n 
var=2.0 
call gnoise(gk,var,k) 
yin(i)=gk 
zkar(i,3)=zkar(i,2) 
zkar(i,2)=zkar(i,1) 
zkma(i,3)=zkma(i,2) 
zkma(i,2)=zkma(i,l) 
zkma(i,1)-yin(i) 
uk(i)=zkma(i,1)*b(i,1)-zkma(i,2)*b(i,2)+zkma(i,3)*b(i,3) 
Zzkar(i,1)=zkar(i,2)*a(i,2)-zkar(i,3)*a(i1,3)+uk(1) 
print *,'--fil-->',afil,theta,r 
print *,'--fil-->',zkar 
yout (i)=zkar(i,1) 
var=save 
return 
end 


$$SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 


subroutine file(yk,theta,step,ic,n,k) 
dimension output (8000) ,theta (10) 

if (k.ne.1) goto 10 

call gendata (output , theta, step, ic,n,k) 
continue 

yk=output (k) 

return 

end 
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$%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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subroutine gendata(output, theta, step,ic,nn,k) 
DIMENSION YSAMP (65536) ,data(4,8000) ,output (8000) 
dimension theta (10) 

OPEN (UNIT=9 , FILE= ' bpsk.dat' ,STATUS='NEW' ) 
n=7000 

amag=1.414 

tcarr=4.8 

tchip=7.0 

ichip=0 

tdata=n 

tdelay=0.0 

1=9999 

tchip=200.0 

do 10 i=1,nn 

tcarr=360.0/theta(i) 

call bpsk(n,amag,tcarr,tchip,ichip,tdata,tdelay,l,ysamp) 
do 20 j=1,n 

if (ic.eq.0) goto 80 

if (i.eq.2) goto 50 

data(i,))=ysamp()) 

goto 60 

if (j.gt.3000) ysamp(j)20.0 
data(i,j)=ysamp(j) 

continue 

continue 

print *,i,'-->',theta(i),tcarr 

continue 

tcarr=360.0/90.0 

tchip=8.0 

call bpsk(n,amag,tcarr,tchip,ichip,tdata,tdelay,1,ysamp) 
do 30 i=1,n 

x=0.0 

do 40 j=1,nn 

x=x+data(j, 1) 

ME (i) =x+yeamp (i) 

write(9,*) x 

continue 

if (ic.eq.0) return 

nl-n-3000 

tchip=200.0 

thetal=theta(2)+step 

tcarr=360.0/thetal 

call bpsk(nl,amag,tcarr,tchip, ichip,tdata,tdelay,1,ysamp) 
eal 

do 70 i=3000,n 
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output (i)=output (i)+ysamp (3) 
j=j+1 

70 continue 

close(9) 

return 

end 


റ 


this programme was modified to suit the martin 
programme 

subroutine bpsk 

& (n,amag,tcarr,tchip,ichip,tdata,tdelay,11,ysamp) 


0000 


THIS PROGRAM GENERATES SAMPLES OF A DIRECT SEQUENCE BI-PHASE 
HIFT 
KEYED SPREAD SPECTRUM SIGNAL. THE "INFORMATION" BITS AND THE 
SPREADING SEQUENCE USED ARE RANDOMLY GENERATED. PARAMETERS 
EQUIRED 
FOR OPERATION ARE: 
N = NUMBER OF SAMPLES GENERATED 
FOR CONSISTENCY WITH USE BY AN FFT 
ALGORITHM, N SHOULD BE AN INTEGER 
POWER OF 2 -- TYPICALLY 1024 
NOTE: IF N>1024 DIM OF YSAMP MUST CHANGE 


MAG = MAGNITUDE OF CARRIER WAVEFORM 
TSDELAY = NUMBER OF SAMPLE TIMES DELAYED FROM 
t = 0 BEFORE BEGINNING SAMPLES 
AN ARBITRARY VALUE THAT ALLOWS SOME 
FLEXIBILITY IN SAMPLING. SHOULD BE 
ZERO FOR SAMPLING AT t=0. 
TDATA = NUMBER OF SAMPLE TIMES IN ONE DATA BIT 
TCHIP = NUMBER OF SAMPLE TIMES IN ONE BIT OF 
SPREADING CODE 
TCARR = NUMBER OF SAMPLE TIMES IN ONE CYCLE OF 
CARRIER FREQUENCY 
TSAMP = DURATION OF A SAMPLE TIME. IN GENERAL 


TSAMP WILL ALWAYS BE = 1.0, SINCE 
TIME SAMPLED VALUES BECOME STATIC 
WHEN STORED AND CAN BE SCALED LATER. 


THE ALGORITHM USED TO GENERATE THE DS-BPSK SIGNAL RELIES UPON 


റടറററററററററററററററററററററററധയററഗററററ 


HE 

BUILTIN FORTRAN RANDOM NUMBER GENERATOR RAND(L) TO PRODUCE THE 
C INFORMATION BIT STREAM AND THE SPREADING SEQUENCE CODE. AT 
THE TIME 
C OF EACH SAMPLE, THESE TWO BITS ARE MODULO 2 ADDED TO PRODUCE 
A CHIP 
C BIT. THIS BIT IS THEN USED TO SHIFT THE PHASE OF THE CARRIER 
BY 
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C 180 DEGREES IF 1 OR BY O DEGREES IF O. THE RESULTANT VALUE OF 
THE : 
C COSINE FUNCTION IS MULTIPLIED BY THE AMPLITUDE OF THE 


WAVEFORM, THEN 

C STORED IN FILE BT.DAT IN THIS FILE DIRECTORY. THE PROCESS IS 
CONTIN- 

C UED UNTIL N SAMPLES ARE GENERATED. NOTE THAT FIRST LINE OF 
BT.DAT 

C CONTAINS N, DSBPSK, DATE AND TIME AS ID FIELD FOR FILE. 

C 

C W.R.TUCKER 9-MAY-83 

e CONVERTED PROGRAM TO ALLOW FOR M-SEQUENCE CHIP 

C M-SEQUENCE GENERATOR REQUIRES DATA FILE FSR.DAT WITH 

C PARAMETERS FOR FEEDBACK SHIFT REGISTER. 

C W.R.TUCKER 4-AUG-83 

C 

C INITIALIZE -- I = INFO BIT NUMBER, J = SP SEQ CODE BIT NUMBER 
C L IS ARGUMENT FOR RAND(L)--A DIFFERENT SEQUENCE MAY BE 

C GENERATED BY INITIALIZING L WITH DIFFERENT VALUES. 

€ 

C Modified by HHL for installation on ULTRIX ECE VAX, 12/90 

e Writes output in file bpsk.dat 

E 


DIMENSION YSAMP(65536) 

INTEGER I,J,L,INFO,CHIP,IFSR,IFBD(100),IVAL(100) 

REAL TDATA, TCHIP, TCARR, TDELAY , TSAMP, MAG, YARG, YSAMP, PI 
PI = 3.14159265 

TSAMP = 1.0 


I = O 
J = O 
L = O 


WRITE(6,900) 
900 FORMAT(' GENERATION OF DIRECT SEQUENCE BPSK SIGNAL', 


A ' IN FILE bpsk.dat') 

WRITE (6,1000) 
1000 FORMAT(' # SAMPLES TO GENERATE = ',$) 
c READ(5, *)N 


print *,n 
WRITE (6,1001) 


1001 FORMAT(' MAX AMPLITUDE OF SAMPLE VALUES (R)= ',$) 
© READ(5, *) MAG 
mag=amag 


print *,mag 
WRITE(6,1002) 
1002 FORMAT(' # SAMPLES PER CARRIER CYCLE (R)= ',$) 
c READ(5, *) TCARR 
print *,tcarr 
WRITE (6,1003) 
1003 FORMAT(' # SAMPLES PER CHIP BIT (R) = ',$) 
c READ(5, *) TCHIP 
print *,tchip 
WRITE(6,1020) 
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1020 FORMAT ( ' ENTER (0):RANDOM CHIP, OR (1):REPEAT 11-51 1118: 
+1,9) 
c READ(5, *) ICHIP 
print *,ichip 
WRITE (6,1004) 
1004 FORMAT(' # SAMPLES PER INFO BIT (R) = ',$) 
c READ(5, *) TDATA 
print *,tdata 
WRITE(6,1005) 
1005 FORMAT(' # DELAYS BEFORE SAMPLING (R)= ',$) 
c READ(5, *)TDELAY 
print *,tdelay 
WRITE(6,1006) 
1006 FORMAT(' RANDOM NUMBER SEED (I4)---» ',$) 


1=11 
C READ(5, *)L 
print *,1 
C Initialize RAND by calling with input non-zero L. 
C Subsequent calls will be with L - O. 
RANDOM-RAND(L) 
e debug 
G WRITE(6,*) ' SEED AND RANDOM NUMBER RETURNED: *,L,RANDOM 
L=0 
WRITE(6,1010) 
1010 FORMAT(/1X,'------------------ WORKING-------------------- ') 
IF (ICHIP .EQ. 0) GO 0025 
OPEN (UNIT = 1, FILE = 'FSR.DAT', STATUS = CED) 


READ (1,500) IFSR 
READ (1,600) (IFBD(K) , K= 1, IFSR) 
500 FORMAT (I3) 
600 FORMAT (I3) 
CLOSE (UNIT = 1) 
DO 5 K= 1 , IFSR 
IVAL(K) = 1 
5 CONTINUE 
DO 100 K = 1,N 
C CHECK TO SEE IF WE NEED TO GENERATE A NEW DATA BIT 
IF ((K+TDELAY)*TSAMP .LT. I*TDATA*TSAMP) GO TO 10 
C IF SO DO TT HERE 
RANDOM = RAND(L) 
debug 
WRITE(6,*) ' IFLAG AND RANDOM NUMBER RETURNED:',L,RANDOM 
INFO = 0 
IF (RANDOM .GT. 0.5) INFO = 1 
C KEEP TRACK OF WHICH DATA BIT WE ARE ON 
I = I+1 
10 CONTINUE 
C NOW CHECK TO SEE IF WE NEED TO GENERATE À NEW CHIP BIT 
IF ((K+TDELAY)*TSAMP .LT. J*TCHIP*TSAMP) GO TO 20 
IF (ICHIP .NE. 1) GO TO 15 
CHIP = IVAL(IFSR) 
CALL MSEQ(IFSR,IFBD,IVAL) 


ററ 
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GO TO 19 
15 CONTINUE 
Sn BP SO DO IT HERE 
RANDOM = RAND(L) 
CHIP = O 
IF (RANDOM GT. 0.5) CHIP = 1 
C KEEP TRACK OF THE CHIP BIT NUMBER 
19 CONTINUE 
J = J+1 
20 CONTINUE 
C NOW WE DETERMINE THE PHASE SHIFT 
PHASE = FLOAT(MOD(INFO + CHIP,2)) 
C COMPUTE THE ARGUMENT FOR THE COSINE FUNCTION 
YARG = 2.0*PI*(1.0/(TCARR*TSAMP) )*(K+TDELAY) + PI*PHASE 
YSAMP(K)= MAG * COS(YARG) 
100 CONTINUE 
C NOW SAVE THE VALUE 
c ranga OPEN (UNIT=1,FILE='bpsk.dat' STATUS='NEW') 
C CALL DATE(BDATE) 


c CALL TIME (CTIME) 

c-ranga  WRITE(1,125)N , MAG ,TCARR, TCHIP, TDATA 

125 FORMAT (15,5X,F4.1,10H*DSBPSK+M2,F7.2,1Hf,F7.2,1Hc,F7.1, 
* ൮99) 


c-ranga WRITE (1,150) (YSAMP(I),I=1,N) 


do 121 i=1,n 
E write(1,*) ysamp(i) 
121 continue 


150 FORMAT(8F16.8) 
CLOSE (UNIT=1) 
return 
END 


SUBROUTINE MSEQ(IFSR,IFDB,IVAL) 
C THIS SUBROUTINE PERFORMS THE SHIFTING OPERATION OF AN IFSR 
STAGE 


C FEEDBACK SHIFT REGISTER, WITH FEED BACK CONNECTIONS AS 
INDICATED 
° PY IFDB. IVAL IS THE INITIAL CONTENTS OF THE FSR AND WILL 
CONTAIN 
E THE FINAL CONTENTS AFTER SHIFTING. 
a W.R. TUCKER 4 AUG 83 

INTEGER IFDB(IFSR), IVAL(IFSR) 

ISUM = O 


DO 10 I = 1, IFSR 

ISUM = IFDB(I) * IVAL(I) + ISUM 
10 CONTINUE 

IBIT = MOD(ISUM,2) 

DO 20 I = 1, IFSR - 1 

IVAL(IFSR + 1 - I) = IVAL(IFSR - I) 
20 CONTINUE 
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IVAL(1) = IBIT 
RETURN 
END 


I2222 2222222222 2222222222222 222222220 
EL IL22I22 222222222222 22222 222222 22222 
22222222 222222222222 2222222222 222222 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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